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Abstract 

We study transition matrices for projected dynamics in the energy-magnetization space, mag- 
netization space and energy space. Several single spin flip dynamics are considered such as the 
Glauber and Metropolis canonical ensemble dynamics and the Metropolis dynamics for three mul- 
ticanonical ensembles: the flat energy- magnetization histogram, the flat energy histogram and the 
flat magnetization histogram. From the numerical diagonalization of the matrices for the projected 
dynamics we obtain the sub-dominant eigenvalue and the largest relaxation times for systems of 
varying size. Although, the projected dynamics is an approximation to the full state space dynam- 
ics comparison with some available results, obtained by other authors, shows that projection in 
the magnetization space is a reasonably accurate method to study the scaling of relaxation times 
with system size. The transition matrices for arbitrary single-spin flip dynamics are obtained from 
a single Monte-Carlo estimate of the infinite temperature transition-matrix, for each system size, 
which makes the method an efficient tool to evaluate the relative performance of any arbitrary local 
spin- flip dynamics. We also present new results for appropriately defined average tunnelling times 
of magnetization and compute their finite-size scaling exponents that we compare with results of 
energy tunnelling exponents available for the flat energy histogram multicanonical ensemble. 

PACS numbers: 02.50.-r,64.60.Ht,05.10.Ln,75.40.Mg,05.50.+q, 02.70.Tt 



I. INTRODUCTION 



The dynamical critical behavior of statistical physics models is a problem that attracts 
considerable attention , y, Q] • From a fundamental point of view one is interested in the 
identification and characterization of the different dynamical universality classes, known to 
be more restricted than the static ones. Different algorithms for canonical ensemble simu- 
lations have been proposed belonging to different universality classes^, 6]. Still, increasing 
relaxation times with system size are a major limitation to the statistical precision of the 
numerical estimates obtained in the simulations. New algorithms^ aiming to estimate the 
number of states of a gi ve» energy, have a,o bee n proposedQ £ Q. These dgorlth™ 
simulate a multicanonical ensemble with the advantage that a single simulation provides 
information on the properties of the system in a wide temperature range. However, such 
algorithms also suffer from slowing down with increasing system size and the study of their 
dynamical properties with simple and efficient methods is essential to ascertain their relative 
performance. 

Many numerical methods have been used to study stochastic dynamics of statistical 
physics models. These methods measure the largest relaxation time of the dynamics, a time 
which increases with system size according to dynamic finite-size scaling theory. The exact 
diagonalization of the transition matrix in the full state space can be done only for very small 
systems. To overcome this limitation, one can instead estimate by Monte-Carlo methods 
the auto-correlation function of the slowest observable in the system, that whose long time 
behavior gives the largest relaxation time. Although this method is free of systematic errors, 
one needs to consider very long simulation runs to get a reasonably small statistical error in 
the auto-correlation function. Several other methods have been used, including a variational 
technique 3|, |4| allowing the estimation of the sub-dominant eigenvalue of the full state-space 
transition matrix. 

Projected dynamics was proposed to study metastability and nucleation in the Ising 
model[10, H, 12, Q, 14, 15]. The idea behind this method is to derive a dynamics in a 



restricted space of one or several variables. Choosing appropriately such variables and ne- 
glecting non-Markovian memory terms one hopes that the resulting approximated Marko- 
vian dynamics is a good approximation to the full state space dynamics. The usefulness 
of the method has been proved in the context of the study of metastability in the Ising 
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model where the direct dynamic Monte-Carlo simulation is unable to bare with the large 



time-scale of the problem 



14J ] . The non-lumpability of the the full state-space transition rate 



matrix with respect to energy and magnetization classification of the states leads to the 

nn 

upcoming of memory terms when projecting the dynamics in these restricted spaces|14j. I16|. 
To recover the Markovian character of the dynamics, these memory terms are neglected and 
the resulting projected dynamics becomes only approximated. 

In this article we study the projected dynamics behavior for the square lattice, nearest- 
neighbor, Ising model, in the energy and magnetization spaces for two local spin flip al- 
gorithms. Namely, the Glauber and the Metropolis et al. 19] critical canonical ensemble 
dynamics, and three multicanonical algorithms: the flat energy-magnetization histogram, 
the flat energy histogram and the flat magnetization histogram dynamics. Although the 
dynamics associated with the transition rate matrices in these restricted spaces are only 
approximate, we show, by comparison with full state space results, that they can be used 
to get reasonably accurate estimates of the dynamical properties. From the numerical di- 
agonalization of these matrices, and the determination of their sub-dominant eigenvalue, 
we compute the largest relaxation times for systems of varying size. The method proposed 
can be applied to other models and other dynamics thus leading to a simple and efficient 
estimation of the scaling with system size of the largest relaxation time. Such studies are 
needed to assess the relative performance of Monte-Carlo simulation algorithms. 

Projected dynamics transition rate matrices were also considered in the context of the 



transition matrix Monte-Carlo 



13, 



181 ]. Using an acceptance probability written in terms 



of the infinite temperature energy space transition matrix it is possible to perform simula- 
tions that visit with equal probability the spectra of energies of the model thus doing flat 
energy histogram simulations. For the case of the Ising model that we consider in this work 
this algorithm is easily generalized to simulations with a flat energy and magnetization his- 
togram. We use this flat energy-magnetization histogram ensemble to numerically estimate 
the infinite temperature transition rate matrix in the space of energy and magnetization 
from which all the results presented in this work were derived. 

For multicanonical algorithms, average tunnelling times between the grou nd-state and 
states with higher energy (for example zero energy) have been considered 20j . It has been 
shown that these tunnelling times may scale differently with system size when we consider 
going up (from a low energy to a high energy) or going down in the energy 21]] . We present 
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new results for average tunnelling times of magnetization, in several multicanonical ensem- 
bles, using projected dynamics, that show a similar behavior and that can be compared with 
results of other authors for tunnelling times in the energy space. 

The new method proposed by us, to study approximately the local dynamics, is efficient 
because: (1) the dynamic exponents estimates are reasonably accurate when compared 
with corresponding quantities obtained by other methods; (2) any, arbitrary, single-spin 
flip dynamics can be studied from a single Monte-Carlo estimation of an infinite temper- 
ature transition matrix in the energy-magnetization space (corresponding to acceptance of 
all the proposed configurations); the consideration of a specific dynamics comes only from 
the weighting of this matrix with the corresponding acceptance probability; (3) the dimen- 
sional reduction achieved by the projection allows the application of matrix diagonalization 
techniques for bigger system sizes. 

The outline of the paper is as follows: In section[TT]we discuss the projection procedure, in 
section ITTT1 we show how the infinite temperature transition matrix is computed from Monte- 
Carlo simulations for different system sizes and we define the projected transition matrices 
for the different ensembles and dynamics considered, in section IIVI we present results for 
the largest relaxation times and the corresponding dynamical exponents, in section |V] we 
define and compute tunnelling times in the magnetization space and their finite-size scaling 
exponents and, finally, in section IVII we summarize our main conclusions. 

II. PROJECTED DYNAMICS 

The Markov chain master equation in the full state space is: 



where a denotes a state of the system, P(a, t) is the probability for the system to be in a 
given state at time t and W(a, a') is the transition rate from state a' to a. In the case of 
an Ising model er = (ax,..., &n) specifies the state of each of iV spins of the system, a i: that 
can take two values, a,i = ±1. The transition rate obeys detailed balance 



dP(a,t) 
dt 



B')P(a', t) - P(a, t)W(a', a)} 



(1) 




(2) 
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relatively to a stationary distribution, P st {a), which we consider to be an arbitrary function 
P st (E(a), M(<?)), of the energy E(a) = j) a i a j (where the sum is over all neighbor 

pairs, and the magnetization M(a) = 

The detailed balance equation can be summed up relative to all a states with a given 
energy E = E(a) and magnetization M = M(a), and all a' states with energy E' = E(a') 
and magnetization M' = M(cr'), to obtain: 



E Pst(&)W (a' , v)5e,e(S)5e>,e(3')5m,m{Z)Smi,m(v') = (3) 

a, a' 

^ Pst(<?')W((T , <t')5e,E(3) &E',E{3>) $M,M(3)$M',M{ff') 

being S a ^ the Kronecker delta. Since the stationary distribution is assumed to be a function 
of the energy and magnetization only, it can be taken out of the summation, giving 

p(E, M)T{E', M'\ E, M) = p(E', M')T{E, M; E' , M') (4) 

where p(E,M) = P st (E , M)Q(E , M) is the stationary probability for the system to have 
energy E and magnetization M, and Q(E, M) is the number of states with energy E and 
magnetization M. In this expression, we have defined, 

T(E', M'\ E, M) = Q ^ M) W (°'> p)$E,E{3)5E> m 5M,Mtff)8 M >,M{ff>) (5) 

as the transition matrix between energy and magnetization states (E,M) and (E',M r ). 

Summing up the master equation in the same way we would obtain the evolution equa- 
tion for the time dependent probability p(E, M, t) for the system to have energy E and 
magnetization M at time t: 

dP{E d ^ t] = E M ' M '> , M '> f ) ~ P( E > M > *)W M '; E > M > 01, (6) 

E',M' 

with a time-dependent transition matrix: 

T(E', M'- E, M- 1) = l M J2 P(a, t)W{a', B)5 e ,e{3)5e>,e{9>) 6 m,m{9) 5 m',m{<}') ( 7 ) 

This time dependent matrix approaches the transition rate matrix in Eq. for large times 
when P(a, i)/p(E, M, t) — > 1/Q(E, M). The so-called projected dynamics neglects this time 



dependence and considers instead the Markov process associated with T(E', M'; E, M): 

dp{E d ^ t] = E [ T ( E > M - E '> M 'M E '> M '> *) - M > l ) T (E\ M'\ E, M)\ . (8) 

E',M' 

Starting with the projection operator technique, in a discrete time formulation, the approx- 
imation can be regarded as equivalent to dropping out some memory terms |X2j] . Note that 
the dynamics of the Markovian process associated with these transition matrices would be 
equivalent to the full state space dynamics if it were lumpablejlfl] with respect to a classi- 
fication of the states in terms of energy and magnetization. However, this is known not to 



be the case for canonical ensemble dynamics 14j] , although the flat magnetization ensemble 



that we study later is lumpable with respect to a magnetization classification of the states. 

Further projection on the energy space can be done by summing for all M and M' the 
detailed balance condition in the E, M space (EqHJ) : 

m E P ^f^-T(E', M'; E, M) = p(E') P{E ' ( m'\ T (E, M; E', M') (9) 

M,M' M,M' ^ ' 

which is a detailed balance relation p(E)T(E', E) = p(E')T(E, E') in the energy space with 
a projected transition matrix 

T(E'; E)=J2 ^r T ^' M '- E > M )- ( 10 ) 

M,M' ' 

Note that for the ensembles where P s tip) depends just on the energy (and not on the 
magnetization) the previous expression can be simplified to: 

T{E ' ] E) = nfE)Yl fi ( E > M ) T (E', M'\ E,M) = ^J2 W ^'' *) S e,e(*Me(«') (11) 

^ ' M,M' 

with Q(E) = J2 M Q(E, M) is the number of states with energy E. If P s t{&) depends on 
energy and magnetization simultaneously the above simplification can not be done. 

The transition matrix T(E, E') can be used to define a Markov chain dynamics in the 
restricted energy space: 

= J2[T(E; E')p(E' 7 1) - p(E, t)T(E'\ E)]. (12) 

E> 

In the same way we can obtain a detailed balance relation in the magnetization space: 

E ^^ T ^'' M < Ej M ) = p ^ E Pi p(M) ) T ( E > M; m ' } ' (13) 



which is a detailed balance relation p(M)T(M' , M) = p(M')T(M, M') in the magnetization 
space with a projected transition matrix 

T(M'; M) = J2 ^WT T ( E 'i M '' E ' M )- ( 14 ) 

E,E' ^ ' 

The transition matrix T(M, M') can be used to define a Markov chain dynamics in the 
restricted magnetization space: 

= ]T[T(M; M>(M\ t) - p(M, t)T(M'; M)]. (15) 

Nevertheless the approximation assumed in the projected dynamics, the detailed balance 
relations satisfied by the transition matrices defined above assure that the long time behavior 
of the related stochastic processeses defined by Eqs. (jSJ), f|T2|) and (|15p are still characterized 
by the correct stationary probability distributions p(E,M), p(E) and p(M), respectively. 

In the following sections, we study single spin flip dynamics in the canonical en- 
semble characterized by the stationary distribution at inverse temperature (3, P s t{&) = 
exp(—/3E(a))/Z as well as three multicanonical ensembles with flat energy-magnetization, 
flat energy and flat magnetization histograms with P s t{&) — 1/Q(E,M), P s t{&) = 1/Q(E) 
and P s t{<j) = l/fl(M), respectively. Note that Q(M) = £l(E, M) is exactly known to 
be Q(M) = ( n+m ) and that an efficient numerical scheme (not used by us in the present 
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work) developed by Beale 22|] allows to compute exactly Q(E) for the two-dimensional Ising 
model for moderate system sizes N. 



III. NUMERICAL CALCULATION OF TRANSITION MATRICES 

We now explain our method to compute numerically the transition matrices 
T(E',M';E,M), T(E',E) and T(M',M) defined in Eqs. ©, $TD]) and ([HD, respectively. 
We start by recalling that for single spin flip dynamics the transition rate W(a', a) can be 
separated in a proposal step and an acceptance step. In the proposal step we choose, with 
equal probability, one of the spins of the system and propose to flip it. Thus a given system 
state may have a non-zero transition rate to iV other system states that differ in the state of 
a single spin. In the acceptance step we accept the proposed configuration with a probability 
a(E', M'; E, M) that we assume depends only on the energy and magnetization of the initial 
and final configurations. 
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Consider the detailed balance relation when we accept all the proposed configurations. 
This is the case, for example, in the Metropolis et al. algorithm at infinite temperature. 
The probability to measure an energy E and magnetization M is then equal to Q(E, M)/2 N 
since all states have equal probability. Thus we can write the relation, 

Q(E, M)T 00 (E', M'\ E, M) = Q(E', M'JT^E, M; E', M') (16) 



231 ] . For a general single spin flip algorithm charac- 



known as the broad histogram equation 
terized by a(E', M'; E, M)we can write, 

T{E', M'\ E, M) = T^E', M'\ E, M)a(E', M'\ E, M). (17) 

The numerical determination of T 00 (£? / , M'; E, M) can be done from the estimator: 

Too{E', M'\ E, M) = NHsm{E ^ M) E N &<> ^wMsi) ( 18 ) 

where the summation is done over the N m configurations generated by the Monte-Carlo 
procedure and N{fik, AE,AM) is the number of configurations with energy E' = E + AE 
and magnetization M' = M + AM that can be obtained from configuration by flipping a 
single spin and the quantity H sim (E, M) is the energy and magnetization histogram of the 
simulation. The estimator in Eq. ffTSj) is related to the average of (N(a, AE, AM))e,m i n 
the constant energy and magnetization ensemble, 

(N(a , AE, AM)) E>M = ^ — g N(a, AE, AM)5 EM3) 5 MM{3) (19) 

= (N{a, AE, AM)8 EMff) 5 MiM{ff) /{P sim {a)n{E, M))) sim 

where P S im(&) is the probability to visit a particular state in the simulation ensemble whose 
averages are denoted by, (. ..) s » m . Since H sim (E,M) = N m P sim (a)Q(E , M) with P sim {v) 
dependent only on E and M, we see that Eq. (TT8]) provides the correct estimator. For the 
two-dimensional square lattice, nearest-neighbor, Ising model each spin can have between 
zero and and four nearest neighbors in the same state of the spin. When this spin flips there 
are five possible energy changes, AE and two magnetization changes, AM. Thus, one needs 
to count the number of spin flips that lead to a energy and magnetization change in each of 
these possible ten classes. 

In this work we have estimated T 00 (£" ', M' '; E , M) by doing transition matrix Monte- 
Carlo simulations in a two-dimensional, nearest neighbor, Ising model of size N = L 2 with 

8 



an acceptance probability given by, a(E' , M'\ E, M) = min ^1, Y^^^j^wWjJ ■ From eqs. 
OH) and f|T7|) we can see that this choice leads to a flat energy and magnetization histogram. 
The algorithm starts with an initial estimate of T oc (i? / , M'; E, M) that is improved as more 
configurations are generated. We have used the n-fold way simulation algorithm of Kalos 



and Lebowitz 



18 



24j an d the number of simulated spin flips per number of spins was 10 8 
for each of the systems studied, L = 3, ...,21,30. Note that, when one considers a n-fold 
way simulation, the histogram of energy and magnetization, H sim (E, M) is the average 
time spent in a given value of energy and magnetization that may differ from the average 
number of hits to a particular energy and magnetization value. In this case, the expression 
(fl8l) should be modified to weight each of the generated configurations with the estimated 
average time spent in these configurations (a small but systematic error arises in the results 
if this weighting is not done). 

The projected transition matrices in the energy-magnetization space, T(E', M'; E, M) 
are obtained from the simulation estimates of T 00 (E' , M' '; E , M) by using eq. (|17|) . We 
consider the following dynamics: (1) the Metropolis canonical ensemble dynamics with 
a(E f , M' '; E , M) = min (1, exp(— (3{E' — E))) (2) the Glauber canonical ensemble dynam- 
ics with a(E',M';E,M) = § (l - tanh(f (E' - E))) (3) the flat energy and magnetiza- 
tion histogram Metropolis dynamics with a(E f , M'; E, M) = min ^1, Y^^J^WId)) > (^) ^ ne 
Metropolis flat energy dynamics also known as entropic sampling with a(E', M'; E, M) = 
min ^1, q ^ j and (5) the Metropolis flat magnetization dynamics with a(E', M'; E, M) = 

min (i'rpro)- 

For the energy magnetization-space with dimension (N + 1) 2 x (N + 1) 2 we have obtained 
results from the diagonalization of T(E', M'; E, M) up to iV = 8 2 . For all the system sizes 
studied we have found the stationary probablities p(E, M) after solving numerically for the 
steady state regime of the system of equations Eq.([Hl)|3j| : 

[T( E , M ; E ', M')p(E', M') - p{E, M)T(E', M'; E, M)] = 0. (20) 

E',M' 

For the flat energy histogram dynamics we need to know Q(E) to construct the corre- 
sponding acceptance probability. This quantity can be obtained from Q(E, M) after the 
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solution of the homogeneous linear system of equations 

£ [Too{E, M; E', M')Sl{E\ M') - Q(E, M)T 00 (E', M'\ E, M)\ = 0. (21) 

E',M' 

Note that it is possible to compute, 

T^E'; E)=J2 ^WeT T ^ E '^ M " E ' M )' ( 22 ) 

M',M \ ' 

and write a(E', M'\ E, M) = min ^1, P^^t^j J for a flat energy histogram ensemble which 
is completely equivalent to a(E', M'\ E, M) = min ^1, 



IV. LARGEST RELAXATION TIMES 

We have considered a discrete time transition matrix defined as j(E, M; E'M') = 
T(E, M; E'M') for (E, M) ^ (E', M') and j(E, M; E, M) = 1 - J2e>,m> T ( E 'i M '' E ' M ) 
for (E,M) = (E',M'). This corresponds to the Markov chain equation p(E,M,t + 1) = 
Y2e' m> l(E, M; E'M')p(E' , M' ,t). The stationary probability distribution corresponds to 
an eigenvector with the largest eigenvalue 1. The second largest eigenvalue, A, determines 
the largest relaxation time in the system, r = ^~ The division by iV is needed in order for 
t to be expressed in units of numbers of Monte-Carlo steps per total number of spins. The 
relaxation times increase as the system size increases as r ~ L z thus being characterized by 
a dynamic exponent, z. 

We studied the critical behavior of the projected dynamics at the critical point of the 
square lattice Ising model, f3J = |ln(l + V2) for Glauber and Metropolis et al. acceptance 
probabilities. Our eigenvalue results for the matrix T(M, M'), \ T ( M > M '} and for the matrix 
T(E, M; E'M'), \T(e,M;E'm') for the Glauber dynamics can be seen in Table [D together with 
the results for the full state space dynamics, X w , obtained from for the Glauber 

dynamics using a variational method. For small systems L = 3,4, 5 we have computed 
yT(E,M;E'M') f rom an e xact enumeration of all the system states and the results are in close 
agreement with the ones obtained from the Monte-Carlo estimation of T^lE, M; E' , M'). 
The eigenvalues, for a given system side, L are close to each other and are observed to obey 
the inequality X w > \nE,M;E>M>) > a t(m ; m')_ 

In Fig. [T] (a) we plot the logarithm of the relaxation time In function of InL, 

for the Glauber dynamics, and the Metropolis et al. dynamics, obtained from the sub- 
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TABLE I: Sub-dominant eigenvalues of transition matrices for different system sides, L, and 
Glauber dynamics. The second column lists values for the matrix W(a,a') taken from Ref. 3j, |4]]. 
The third and fourth columns are our results for the matrices T(M; M') and T(E, M; E'M') 
respectively. 



L 


X w Rei [3, £ 


A T(M;M') 


\T(E,M;E',M') 


3 


0.997409385126011" 


0.9973901755 


0.99740630184576° 








0.9974063007 


4 


0.999245567376453° 


0.9992429803 


0.99924409354918° 








0.9992441209 


5 


0.999708953624452° 


0.9997066202 


0.99970673172786° 








0.9997067351 


6 


0.9998657194 


0.9998635780 


0.9998637800 


7 


0.9999299708 


0.9999281870 


0.9999284453 


8 


0.9999600854 


0.9999586566 


0.9999589090 


9 


0.9999756630 


0.9999744986 




10 


0.9999843577 


0.9999834244 




11 


0.9999895056 


0.9999887396 




12 


0.9999927107 


0.9999921039 




13 


0.9999947840 


0.9999942741 




14 


0.9999961736 


0.9999957520 




15 


0.9999971315 


0.9999967823 




16 


0.9999978080 


0.9999975119 




17 


0.9999982987 


0.9999980505 




18 


0.9999986606 


0.9999984474 




19 


0.9999989315 


0.9999987550 




20 


0.9999991370 


0.9999989750 




21 


0.9999992955 


0.9999991723 




30 




0.9999998016 





°Exact 
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dominant eigenvalue of T(M;M'), together with the full state space results of [3j, |4j and 
also te for the Glauber dynamics, obtained from T(E;E'). The fitted straight lines were 
obtained neglecting data for L < 15 and have slopes z GL = 2.02, z Met ' = 2.00 and z = 2.18. 
To estimate a reliable value of the exponent z a careful analysis is needed taking into 
account corrections to scaling. In [Jjthe authors report their best estimate z = 2.1660(10) 
excluding the Domany conjecture z = 2 with a logarithmic factor r ~ L 2 (l + 51nL) 25]. 
Further analysis 2q| . by other authors, of the data of Nightingale and Blothe was not able 
to categorically exclude the validity of the Domany conjecture. The logarithmic dependence 
of te on system size, in accordance with previously reported results 17J] , can be seen in Fig. 
CD (a). We have not tried to make a detailed analysis of our results in order to have precise 
estimates of the dynamic exponent from the magnetization projected dynamics. However, 
in Fig. [T] (b) we show the local slope z = i^t^^^tj^ as a function of L~ 2 which is the first 
order finite size correction to the leading behavior considered in 3j, |4|, t ~ L z (l + bL~ 2 ). 
The results of the extrapolation to the infinite system size limit shown in figure Fig[T] (b) are 
z gi. _ i_gg ; z Met. _ 2.00 and z = 2.165. The results for z GL and z Met ' seem to be consistent 
with z GL ~ z Met - ~ 2. 

In Fig. Ufa) we plot the relaxation times obtained from T(M; M') for the Metropolis 
et al. dynamics in the flat energy-magnetization ensembles, flat energy ensemble and flat 
magnetization ensemble. The fitted straight lines have slopes given by z^m = 2.11, z 1 ^ = 
2.69 and z^ = 1.99, respectively. In Fig. [2(b) we plot the relaxation times obtained from 
T(E; E') for the Metropolis et al. dynamics in the same ensembles. The slopes of the fitted 
straight lines are z% M = 2.14, z% = 2.13 and zf { = 1.99, respectively. We neglected the 
data for L < 15 in all the fits shown in Figs. [2] (a) and (b). 

In Fig. [S](a) we show the estimation of the dynamic exponent from the local slopes of 
the relaxation time plots in the magnetization projected space (Fig. a)) as a function of 
L~ 2 . The estimates for z^ (flat energy histogram) in Fig. [3](a) seem to have converged for 
the system sizes studied. The average value for system sides L > 15 is z^ = 2.68(2). This 
result is compatible with the available [2 71] result z = 2.80(13) obtained by a Monte Carlo 
estimate of the convergence time of the time-dependent energy histogram to the stationary 
flat distribution of the energy. The estimates for Zj% (flat magnetization histogram) plotted 
in the same figure show a dependence with system size approaching a value close to 2 for 
infinite system sides. The infinite size extrapolation for z^m energy and magnetization 
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L~ 2 



G. 1: (a) Values of lnr^ for Glauber dynamics obtained from the results reported in Ref. 
J] (o) and values obtained from T(M;M') (+) as a function of InL. We also plot In tm for 
the Metropolis et aJ. dynamics (□) and r^; for the Glauber dynamics obtained from the matrix 
T{E;E') (o). The fitted straight lines were obtained neglecting data for L < 15 and have slopes 



Ml. 



2.02, z 



Met. 



2.00 and z = 2.18. (b) Estimations of the dynamic critical exponent from the 



local slopes of the graph in (a) as a function of L 2 . The symbols are as in (a). The extrapolated 
exponents are z GL = 1.991, z Met = 2.00 and z = 2.165. 

histogram) gives a value slightly larger than 2, z^ M = 2.08. For this multicanonical dynam- 
ics the results show an even-odd effect and it is important to do separate estimates for even 
and odd system sides. In Fig. [3](b) we also show the estimation of the dynamic exponent 
from the local slopes of the plots (Fig. Mjo)) for the energy projected dynamics. The infinite 
size extrapolation for z§ M for even and odd system side are very close to each other and 
given by, 

z e m — 2.07. The extrapolations for for odd system side and even system side 
give Ze = 2.07 and 1.99, respectively. The difference between these two estimates may be a 
sign of the presence of corrections to scaling not properly accounted by our analysis. The 
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FIG. 2: Relaxation times obtained for the Metropolis et al. dynamics in the flat energy- 
magnetization ensemble (o), flat magnetization ensemble (*) and flat energy ensemble (x). In 
(a) we plot tm obtained from T(M;M') and the fitted straight lines have slopes z^f M = 2.11, 
z¥ = 2.69 and zH = 1.99. In (b) we plot te obtained from T(E;E') and the fitted straight lines 
have slopes z^ M = 2.14, z§ = 2.13 and z?, = 1.99. Both in (a) and (b) data for L < 15 were 
neglected in the fits. 

estimates for zfj show a size dependence that is compatible with a value close to 2. 

The full state transition rate matrix W, in the flat magnetization ensemble is lumpable 
with respect to the classification of the states according to their magnetization. Conse- 
quently, the result zj% = 2.00 does not suffer from the approximation inherent to the pro- 



jection procedure. A sufficient and necessary condition for lumpability[16|, is that the total 
probability to go from a state belonging to a given magnetization class to another class with 
different magnetization is the same for every state in the starting class. For each state in the 
starting class with magnetization M there are n± states in the final class M ± 2 where n± is 
the number of up/down spins in the initial configuration. The probability to move to each 
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L- 2 L-2 

FIG. 3: (a) Dynamic exponent estimates z^l M {o), z|f(x), from local slopes of the plots 

shown in Fig. [2] (a) as a function of L~ 2 . For the estimates for the energy-magnetization flat multi- 
canonical ensemble we made separate estimations for even and odd L system sides. The infinite sys- 
tem size extrapolations give, z^ M = 2.08, = 2.68 and zjfa = 2.00. (b) z§ M (o), z§(x), zf f (*) 
from local slopes of the plots shown in Fig. [2] (b) as a function of L~ 2 . For the estimates for the 
flat (E, M) ensemble and for the flat energy ensemble we made separate estimations for even and 
odd L system sides. The infinite system size extrapolations give, z§ M = 2.07, z§ = 2.07 and 1.99 
for odd and even L, respectively, and zfj = 2.00. 

of these final states in the final class has a constant value that depends only on the initial M 
and on the final M±2. All the states in the starting class have the same number of up spins 
and down spins so the probability to move to M±2 is the same for every state in the starting 
class. The matrix T(M; M') is a tridiagonal symmetric matrix with matrix elements given 
by, T(M + 2;M) = T(M;M + 2) = ^±±i for M < 0, andT(M + 2;M) = T(M;M + 2) = ^ 
for M > 0. 
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V. MAGNETIZATION TUNNELLING TIMES 



As a measure of performance for multicanonical methods the average tunnelling times 
were introduced [20]. These tunnelling times measure the time required to sample the whole 
phase space and scale with system size differently than the relaxation time. It was shown 
that it is important to distinguish between tunnelling from ground-sate to maximum energy, 



2l| 



the up direction, and from the high energy to the ground-state, the down direction 

All the tunnelling times reported by us are calculated for the projected dynamics as- 
sociated with T(M; M'). We calculate the average time, r t for the system to go from 
magnetization M = —N to M = N. We also consider two other average times: The time r u 
for the system to go from M = —N to zero magnetization, and the time Td for the system 
to go either to M — +N or M = —N when it starts from M — 0. This definition of t u and 
Td apply only to systems with even L (and N) such that M = is an accessible value of the 
magnetization. 

The tunnelling times above defined obey the relation t u + tj = r t /2 that follows from the 
following simple argument: For the system to go from M = —N to M = N it has to reach 
M = at some point. It will do so for the first time using an average time r u . Then with 
probability 1/2 it will reach for the first time M = N and the tunnelling time would be r u +Td 
or it will return to M — —N and it will reach later M = N taking a time r t . Consequently 
the tunnelling times obey the relation \{j u + r d ) + \r t = r t . This argument uses the fact 
that the matrix T(M; M') has the symmetry property T(M±2; M) = T(-M=f2; -M) and 
so the walk along positive values of the magnetization has the same statistical properties of 
the walk along negative values of the magnetization. 

The time to go from M = —N to M = N can be easily computed taking advantage of 
the fact that T(M; M') is non zero only when M = M' ± 2 and M' = M. If we do not allow 
transitions from M = N to M = N — 2 the M = N becomes an absorbing site for every 
walk along the magnetization axis meaning that it will end there upon a first visit. Defining 
h(M) as the average time spent at magnetization value M [jjj], we can write: 

h(M - 2)T(M; M - 2) - h(M)T(M - 2; M) = 1 (23) 

which means that the difference between the average number of jumps in the positive direc- 
tion (M — 2 — > M) and the average number of jumps in the negative direction (M — ► M — 2) 
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should be equal to one since the system will eventually reach M = N by moving one time 
in excess in the positive direction through the bond connecting the sites M — 2 and M. At 
M = N there are no jumps in the negative direction and so h(N — 2)T(N; N — 2) = 1. It 
is then simple to calculate h(M) and the average tunnelling time for the system to go from 
M = -N to M = N is given by,r 4 = Em=-/ H m )- 

The time t u to reach for the first time M = starting from M = —N is obtained using 
the recursion ([21?]) together with the equation h(— 2)T(0; —2) = 1 to get r u = ^m=-at h(M). 
Finally, the average time required to start from M = and reach for the first time either 
M = — iV or M = N, is obtained from the recursion 

h(M)T(M - 2; M) - h(M - 2)T(M; M — 2) = 1 (24) 

with a modified rate T(-2; 0) equal to T(-2; 0) + T(2; 0) and h(-N+2)T(-N; -N + 2) = 1. 
The average time r d is then given by = Sm=-at+2 h(M). The average tunnelling times 
obtained by this method could also have been obtained from the calculation of the probability 
of first visit to the absorbing site that can be computed from the eigenstates and eigenvectors 



of the associated absorbing Markov chain matrix ( see [2l|). 

The tunnelling times are characterized by dynamic exponents[20], r t ~ L d+Zt , r u ~ L d+Zu , 
Td ~ L d+Zd . The relation between these tunnelling times imply that Zt is equal to the biggest 
of the two exponents, z u and Zd, z t = max(zd, z u ) . Note that the tunnelling times reported 
by us are measured in units of lattice sweeps and not in units of single site updates. 

In Fig. H](a) we show the size dependence of the tunnelling times, r t , r u and for the 
Metropolis et al. dynamics in a flat magnetization-energy histogram ensemble obtained 
from the matrix T(M; M'). We see that r t ~ ^> t u . The scaling exponent of r u obtained 
from the plot is z u = 0.15 which predicts a scaling t u ~ L 215 close to the exponent of the 
relaxation time z^ M = 2.11 reported in the previous section. This behavior is similar to the 



one found in j2l| where t u (in the energy space) was found to scale like the relaxation time 
of the system. For the other two exponents we have obtained z t ~ Zd ~ 2.12. Note that for 
a random walk in the magnetization axis a value for these exponents equal to is expected. 
In Fig. IU(b) we show the local slopes for the plots in Fig. Ufa) as a function of L~ 2 . The 
estimates for z t and Zd seem to follow a straight line predicting an infinite system value 2.11 
and 2.08, respectively. The infinite size extrapolation for z u is 0.14. 

In Fig. 0(a) we show the size dependence of the tunnelling times, r t , t u and r d for 
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FIG. 4: (a) Tunnelling times, (o), t u (+) and (*) as a function of system size L for the 
Metropolis et al. dynamics in a fiat magnetization-energy histogram ensemble obtained from the 
matrix T(M;M'). The fitted straight lines were obtained neglecting data for L < 15 and have 
slopes, zt = 2.12, z u = 0.15 and = 2.12, respectively. In (b) we plot the corresponding local 
slopes as a function of L~ 2 . Even and odd system sides were treated separately. The infinite 
system extrapolation gives, z% = 2.11, z u = 0.14 and z& = 2.08. 

the Metropolis et al. dynamics in a flat energy histogram ensemble obtained from the 
matrix T(M; M'). The slopes of the fitted straight lines give z t = 0.69, z u = 0.64 and 



Zd = 0.63. The result for z t can be compared with the value 0.78 reported in Ref. [28j and 



the value 0.743(7) reported in Ref. [20] by measuring average times for energy excursions. An 
exponent z u = 0.6, also obtained from Monte-Carlo estimates of energy tunnelling times was 
previously reported [29| in very good agreement with our result. In Fig. E^b) we make infinite 
size extrapolations giving, z t = 0.70 and 0.65 for odd and even system sides, respectively, 
z u = 0.63 and Zd = 0.66. 

Finally, we consider the Metropolis et al. dynamics for the flat magnetization histogram 
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FIG. 5: (a) Tunnelling times, Tt (o), r n (+) and r d (*) as a function of system size L for the 
Metropolis et al. dynamics in a flat energy histogram ensemble obtained from the matrix T(M; M'). 
The fitted straight lines were obtained neglecting data for L < 15 and have slopes, Zt = 0.69, 
z u = 0.64 and z d = 0.63, respectively. In (b) we plot the corresponding local slopes as a function 
of L~ 2 . For the zt estimates even and odd system sides were treated separately. The infinite size 
extrapolations are, Zt = 0.70 and 0.65 for odd and even system sides, respectively, z u = 0.63 and 
z d = 0.66. 

ensemble. For this case it is possible to compute analytically the tunnelling times from the 
recursion relations given above, Eqs. (I23|I24I) and the knowledge of the matrix T(M;M'). 
The analytical results are, t u = N/2, r t = (N + l)H(N/2) and r d = \ ((JV + l)H{N/2) - N) 
where, H(n) = 5^fc=il/^ ^ s ^ ne Harmonic number. Using the known asymptotic result, 
for large n, H(n) ~ (Inn + 7) where 7 = 0.5772156649... is the Euler constant we have 
asymptotic expressions for the tunnelling times that predict, r t /N ~ t^/N ~ lniV and the 
tunnelling exponents are, z t = z d = z u = 0. In Fig. El(a) we compare the numerical results 
for the tunnelling times, r t , r u and r d with the analytical results. Note that, because of the 
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FIG. 6: (a) Tunnelling times, (o), t u (+) and (*) as a function of system size L for the 
Metropolis et ai. dynamics in a flat magnetization histogram ensemble obtained from the matrix 
T(M; M'). The lines are the analytical asymptotic results given in the text. 



logarithmic dependence of Tj and r u the estimates for the exponents Zt and Zd that we could 
obtain for the slopes of the data shown in Fig. E](a) give effective values around 0.35 that 
would slowly approach zero only if larger systems were considered. 

From the three multicanonical ensembles studied we see that the flat magnetization en- 
semble is the one with smaller tunnelling exponents and relaxation time exponent. Recently, 
it was shown that it is possible to optimize the ensemble in multicanonical simulations such 



that the tunnelling exponent, z t is also reduced to zero 



28. 



30|. 



VI. CONCLUDING REMARKS 



We have shown that projected dynamics in the magnetization space is a reasonably good 
approximation to the full state space single spin flip dynamics studied in this work: canonical 
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ensemble Glauber and Metropolis et al. dynamics and three multicanonical ensemble dy- 
namics with flat energy-magnetization, flat energy and flat magnetization histograms. The 
energy projected dynamics is generally a worse approximation being not able to preserve 
the power-law size increase of the relaxation time for canonical ensemble dynamics. From 
all the studied dynamics only the flat energy histogram dynamics show a z exponent clearly 
larger than 2 and near 2.7. For the case of the flat magnetization histogram the projection 
in the magnetization space is exact and it is possible to obtain analytical results for the 
tunnelling times predicting a zero value for the exponents, z t , Zd and z u . The tunnelling 
exponents, Zt (and Zd) for the energy and magnetization flat histogram ensemble are much 
bigger, z t — Zd ~ 2 and larger than the exponent z u ~ 0. For the flat energy histogram 
dynamics these three exponents are not very different and the estimates fall between the 
values, z u ~ 0.63 and z t ~ 0.70 for odd system sides. These results were obtained from the 
tunnelling properties of the projected dynamics in the magnetization space that were found 
to be in rough agreement with ones obtained by independent methods for excursions in the 
energy space for the flat energy multicanonical ensemble. 

Finally, the results show that the evaluation of the relative performance of single-spin 
flip dynamics in Ising like models can be done very efficiently by studying the projected dy- 
namics in the magnetization space: the approximation gives reasonably accurate dynamic 
exponents; any, arbitrary, single-spin flip dynamics can be studied from Monte-Carlo es- 
timations of Too(E, M; M') for several system sizes in the energy-magnetization space 
and the large dimensional reduction achieved by the projection in the magnetization space 
allows the application of matrix diagonalization techniques for bigger system sizes. 

Furthermore, the application of projection methods to cluster dynamics in Ising models 
and also to other models projected along their slowest mode may be of considerable interest. 
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